!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine damping_zone (ro, mx, my, mz, bx, by, bz, en,                                                                           &
   ro0, mx0, my0, mz0, bx0, by0, bz0, en0, x, y, z, nx, ny, nz)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny, nz

   real(8), intent(in) :: x(nx), y(ny), z(nz)
   real(8), dimension(nx, ny, nz), intent(in) :: ro0, mx0, my0, mz0, bx0, by0, bz0, en0
   real(8), dimension(nx, ny, nz), intent(inout) :: ro, mx, my, mz, bx, by, bz, en

   real(8) :: damps_x, dampe_x, dampa_x, dampb_x
   real(8) :: damps_y, dampe_y, dampa_y, dampb_y
   real(8) :: damps_z, dampe_z, dampa_z, dampb_z
   real(8), dimension(nx, ny, nz) :: damp
   real(8), dimension(nx, ny, nz) :: pr, pr0
   
!----------------------------------------------------------------------------------------------------------------------------------|
   damps_x = -0.4d10
   dampe_x = -0.5d10
   dampa_x = 6.d0 / (dampe_x - damps_x)
   dampb_x = -(damps_x + dampe_x) / 2.d0
   damps_y = 0.4d0
   dampe_y = 0.5d0
   dampa_y = 6.d0 / (dampe_y - damps_y)
   dampb_y = -(damps_y + dampe_y) / 2.d0
   damps_z = 0.4d10
   dampe_z = 0.5d10
   dampa_z = 6.d0 / (dampe_z - damps_z)
   dampb_z = -(damps_z + dampe_z) / 2.d0
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx
      damp(i, j, k) = (-tanh (dampa_x * (x(i) + dampb_x)) + 1.d0) / 2.d0 *                                                         &
                      (-tanh (dampa_y * (y(j) + dampb_y)) + 1.d0) / 2.d0 *                                                         &
                      (-tanh (dampa_z * (z(k) + dampb_z)) + 1.d0) / 2.d0
      pr0(i, j, k) = gmm1 * (en0(i, j, k) -                                                                                        &
         (mx0(i, j, k) * mx0(i, j, k) + my0(i, j, k) * my0(i, j, k) + mz0(i, j, k) * mz0(i, j, k)) / ro0(i, j, k) / 2.0d0 -        &
         (bx0(i, j, k) * bx0(i, j, k) + by0(i, j, k) * by0(i, j, k) + bz0(i, j, k) * bz0(i, j, k)) / 2.0d0)
      pr(i, j, k) = gmm1 * (en(i, j, k) -                                                                                          &
         (mx(i, j, k) * mx(i, j, k) + my(i, j, k) * my(i, j, k) + mz(i, j, k) * mz(i, j, k)) / ro(i, j, k) / 2.0d0 -               &
         (bx(i, j, k) * bx(i, j, k) + by(i, j, k) * by(i, j, k) + bz(i, j, k) * bz(i, j, k)) / 2.0d0)
   enddo
   enddo
   enddo

   do k = 1, nz
   do j = 1, ny
   do i = 1, nx
      en(i, j, k) = en(i, j, k) - pr(i, j, k) / gmm1
      ro(i, j, k) = ro0(i, j, k) * (1.d0 - damp(i, j, k)) + ro(i, j, k) * damp(i, j, k)
      mx(i, j, k) = mx0(i, j, k) * (1.d0 - damp(i, j, k)) + mx(i, j, k) * damp(i, j, k)
      my(i, j, k) = my0(i, j, k) * (1.d0 - damp(i, j, k)) + my(i, j, k) * damp(i, j, k)
      mz(i, j, k) = mz0(i, j, k) * (1.d0 - damp(i, j, k)) + mz(i, j, k) * damp(i, j, k)
      pr(i, j, k) = pr0(i, j, k) * (1.d0 - damp(i, j, k)) + pr(i, j, k) * damp(i, j, k)
      en(i, j, k) = en(i, j, k) + pr(i, j, k) / gmm1
   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
end subroutine
